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ABSTRACT 

Despite the importance of their size evolution in understanding the dynamical evolution of globular 
clusters (GCs) of the Milky Way, studies are rare that focus specifically on this issue. Based on the 
advanced, realistic Fokker-Planck (FP) approach, we predict theoretically the initial size distribution 
(SD) of the Galactic GCs along with their initial mass function and radial distribution. Over one 
thousand FP calculations in a wide parameter space have pinpointed the best-fit initial conditions 
for the SD, mass function, and radial distribution. Our best-fit model shows that the initial SD of 
the Galactic GCs is of larger dispersion than today's SD, and that typical projected half-light radius 
of the initial GCs is ~4.6 pc, which is 1.8 times larger than that of the present-day GCs (~2.5 pc). 
Their large size signifies greater susceptibility to the Galactic tides: the total mass of destroyed GCs 
reaches 3-5 xlO 8 M Q , several times larger than the previous estimates. Our result challenges a recent 
view that the Milky Way GCs were born compact on the sub-pc scale, and rather implies that (1) 
the initial GCs are generally larger than the typical size of the present-day GCs, (2) the initially large 
GCs mostly shrink and/or disrupt as a result of the galactic tides, and (3) the initially small GCs 
expand by two-body relaxation, and later shrink by the galactic tides. 

Subject headings: Galaxy: evolution - Galaxy: formation - Galaxy: kinematics and dynamics - globular 
clusters: general - methods: numerical 



1. INTRODUCTION 

Whereas the present-day mass functions (MFs) of glob- 
ular cluster (GC ) systems, which are nearly universal 
amon g galaxies (|Brodie fe StradeH 120061 : Uornan et al.1 
12007ft . are approximately log- normal with a peak mass 
M p « 2xl0 5 M©, the MFs of the young massive 
star cluster (YMC) systems follow a sim ple power-law 
distribution itWhitmore fe Schweizer1 [l995: Zh ang fe Fall 



119991: Ide Griis et al.l 120031 among others). Motivated 
by such a difference between GCs and YMCs, nu- 
merous studies have examined the dynamical evolu- 
tion of the GC MFs to determine whether the initial 
MFs of GC systems resemble those of YMC systems 
(iGnedin fe Ostriker l )1997t iBaumgardtl 119981: iVesperinil 
19981: iFall fe Zhangl (200l[ iParmentier fe Gilmorel l2007t 
Shin, Kim, & Takahashi 2008, among others). In partic- 
ular. iShin. Kim, fe Takahashil 1)20081 " Paper I hereafter) 
surveyed a wide range of parameter space for the initial 
conditions of the Milky Way GCs, and considered virtu- 
ally all internal/external processes: two-body relaxation, 
stellar evolution, binary heating, galactic tidal field, ec- 
centric orbits and disc/bulge shocks. They found that 
the initial GC MF that best fits the observed GC MF 
of the Milky Way is a log-normal function with a peak 
at 4xl0 5 Mq and a dispersion of 0.33, which is quite 
different from the typical MFs of YMCs. 
Using the outcome of Y-body calculations, 
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IGieles fc Baumgardt] (|2008f ) found that the aspect 
of mass loss in GCs varies with the tidal filling ratio 
5ft = rh/rj, where ry t is the half-mass radius and rj 
is the Jacobi radius. More specifically, the mass loss 
of GCs in the "isolated regime" (5ft < 0.05) is driven 
mostly by the two-body relaxation, which induces the 
formation of binaries in the core and causes GCs to 
expand. On the other hand, the mass loss of GCs in the 
"tidal regime" (5ft > 0.05) is influenced by the galactic 
tides as well, which enables stars in the outer envelope 
to easily escape (evaporation). Thus, the cluster size 
(rh) is as important as the cluster mass (M) and the 
galactocentric radius (Rg) in determining the dynamical 
evolution of GCs. 

Can YMCs tell us something about the typical ini- 
tial size of the Milky Way GC system? Observations 
show that the projected half-light radius Rh of YMCs 
(ages up to 100 Myr) in the local group ranges be- 
tween ~2 and ~30 pc with a mean valu e of ~8 pc 
(|Portegies Zwart. McMillan, fc Gielesll20ldl ). which is a 
few times larger than that of the present-day Milky Way 
GCs, Rh ~ 2.5 pc. However, GCs could have formed in 
different environments and/or by different mechanisms 
from the YMCs. 

Perhaps the best way to estimate the typical size of 
the GCs is to trace them back to their initial state by 
calculating their dynamical evolution. In this paper, we 
study the dynamical evolution of the Galactic GCs and 
identify the most probable initial conditions not only for 
the MF and radial distribution (RD), but also the size 
distribution (SD). Using the same numerical method and 
procedure as in Paper I, we perform Fokker-Planck (FP) 
calculations for 1152 different initial conditions (mass, 
half-mass radius, galactocentric radius and orbit eccen- 
tricity), and then search a wide-parameter space for the 
most probable initial distribution models that evolve into 
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Fig. 1. — Distribution of the Galactic "native" (see the text for 
definition) globular clusters in the L-Rq space (a), R^-Rq space 
(b), and Rh-L (c). Data are from the compilation by Harris (1996). 

the present-day Galactic GC distributions. 

The paper is organized as follows. Section 2 describes 
the properties of the observed GCs, again which we com- 
pare our model results. Section 3 presents models and 
initial conditions for FP calculations, and Section 4 an- 
alyzes the aspects of the size evolution of GCs. We syn- 
thesize our FP results in Section 5 to construct the GC 
system, and examine common features of the best-fit MF, 
RD, and SD models in Section 6. We discuss characteris- 
tics of the final best-fit SD models of the Galactic GCs in 
Section 7. Finally, conclusions are presented in Section 
8. 

2. PRESENT-DAY GC PROPERTIES 

When comparing FP calculations to the present-day 
Galactic GCs, we consider the "native" GCs only, 
i.e., "old" halo and bulge/disc clusters, which are be- 
lieved to be created when a protogalaxy collapses while 
"young" halo clusters are thought to be formed in exter- 
nal satellite galax ies (Zinn 1993; iParmentier etai]l2000t 
iMackev & ' van den Berghl |2005| ). Our native GC can- 
didates do not include six objects that belong to the 
Sagittarius dwarf, seven objects whose origins remain un- 
known, two objects that have no size information, and 
fifteen objects that are thought to be the remna nts of 
dwarf galaxies (jLee. Gim. fc Casetti-Dinescul l2007h The 
total number of our present-day Galactic native GCs is 
93, and their observed properties, such as luminosity L, 
Rh , and Rq, we re obtained from the database compiled 
bv lHarrislll996l) . 

Figure [1] shows scatter plots between observed L, Rh, 
and Rq values for the 93 Galactic native GCs. The L, 
Rh, and Rq values range between 3.9x 10 3 -5.0x 10 5 Lq, 
0.3-16 pc, and 0.6-38 kpc, where the mean values are 
located at 7.2xl0 4 L Q , 2.5 pc, and 4.1 kpc, respectively. 
The correlation between Rh and Rq is tighter than the 
other two correlations (see Figure QJ). This tight Rh~ 
Rg correlation could be just a result of the initially tight 
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Fig. 2. — Comparison of evolution between iV-body simula- 
tions and our FP calculations for GCs with initial conditions of 
M = 10 4 Mq, Rq = 8.5 kpc, and rj, = 1, 3, and 5 pc on circu- 
lar orbits. TV-bo dy simulations were performed using Nbody4 code 
(Aarseth 2003), and mass loss by stellar evolution was not consid- 
ered in these test calculations (both JV-body and FP). values of 
the two models agree well within ~ 20% during the entire cluster 
lifetimes. 

correlation between Rh and Rg, or it could be due to the 
preferred disruption of larg e GCs near the Galactic center 
i|Vesperini fc HeggidH997t iBaumgardt fc Makinoll2003[) . 
Another possible cause is the expansion of initially small 

2/3 

GCs up to rj, which is roughly proportional to i? G for 
a given GC mass. One of the goals of this paper is to 
determine which of these possibilities is more feasible. 

Previous studies on the evolution of the GC system as- 
sumed a certain constant mass-to-light (M/L) ratio, and 
converted the observed L to M when comparing their nu- 
merical values with observations. But the conversion of 
GC luminosity function (LF) to GC MF using a constant 
M/L ratio may lead to MFs in error because low-mass 
stars, which have higher M/ L ratios than the high-mass 
stars, preferentially evaporate from the cluster and this 
causes the M / L ratio of the cluster to evolve with time 
(|Kruiissen &: Portegies Z wart 2009[). For the same rea- 
son, there is not a linear relationship between Rh and 
r/j among different GCs. Thus, we transform M to L, 
instead of L to M, using the stellar mass-lumin osity re- 
lation of the Padova model (Marig o et al.l l2008h with a 
metallicity of [Fe/H] = —1.16, which is the mean value 
for the Galactic native GCs. Our FP calculations, which 
will be described later, show that the present-day GCs 
can have M/L ratios ranging between 1.2 and 2.5 and 
Th/Rh ratios ranging between 1.0 and 2.5. 

We use dynamical properties such as M and r/j when 
constructing the initial distributions of the Galactic GC 
system and when calculating the dynamical evolution, 
while observed quantities, L and Rh, are used when com- 
paring our FP results with the observations. 

3. MODELS AND INITIAL CONDITIONS 
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We adopt the anisotropic FP mod el used in Paper 
I, whi ch was originally developed by iTakahashi fc Lee] 
(|2000t and references therein). The model integrates 
the orbit-averaged FP equation of two (energy- angular 
momentum) dimensions and considers multiple stellar 
mass components, three-body and tidal-capture binary 
heating, stellar evolution, tidal fields, disk/bulge shocks, 
dynamic al friction, a nd re alistic (eccentric) cluster or- 
bit (see iKim fc Led ([T999) for the tidal binary heat- 
ing and Paper I for the detailed implementation of dy- 
namical friction and realistic orbits). The model imple- 
ments the Al ternating Direction Implicit (ADI) method 
developed bv lShin fc Kim! ([2007) for integrating the two- 
dimensional FP equation with better numerical stability. 

Parameters for our FP survey are the following four 
initial cluster conditions: M, r^, apocenter distance of 
the cluster orbit R a , and cluster orbit eccentricity e. We 
choose eight M values from 10 3 5 to 10 7 Mq, six values 
from 10 _1 to 10 pc, and six R a values from 10° to 
10 167 kpc, all equally spaced on the logarithmic scale. 
For the eccentricity, we choose e = 0, 0.25, 0.5, and 0.75. 
We perform FP calculations for all possible combinations 
of these four parameters, thus the total number of cluster 
models considered in the present study amounts to 1152. 

For the initial stellar mass function (IMF) within each 
cluster, we adopt the model developed by Kroupa (2001) 
with a mass range of 0.08-15 M Q , which is realized by 15 
discrete mass components in our FP model. Each mass 
co mponent follows the s tellar evolution recipe described 
bv IS dialler et al.l (J1992). The stellar density and veloc- 
ity dispersio n distributio ns within each cluster follow the 
King model (King 1966) with a concentration parameter 
Wq = 7 and with neither initial velocity anisotropy nor 
initial mass segregation. We use only one value for Wo, 
thus the tidal cut-off radius r t of the King profile is pro- 
portional to rh , while r j varies depending on M and Rq . 
Therefore, the Roche lobe filling ratio (r t /rj) and 5ft of 
our FP models are functions of r^, M, and Rg- 

The aspects of mass and size evolution from our FP 
model are in a good agreement with those from A^-body 
methods. A comparison of mass evolution between our 
FP calculations and the T V -body simulations performed 
by iBaumgardt fc Makinol (|2003f ) for clusters on eccen- 
tric orbits with initial masses larger than 10 4 Mq shows 
good agreement of cluster lifetimes within ~25%. For a 
comparison of size evolution, we run a set of A ^-body 
simulations using Nbody4 code (jAarseth 1 120031 ) with 
M = 10 4 M , Rq = 8.5 kpc, and = 1, 3, and 5 pc 
(these correspond to 5ft = 0.04, 0.11, and 0.18), and find 
that the evolutions of the two models agree well within 
~ 20% during the entire cluster lifetimes (see Figure [2]). 

Due to the expulsion of the remnant gas from star 
formation in the pre-gas-expulsion cluster, some of 
the low-mass pre-gas expulsion clusters can quickly 
disrupt, and even the surviving low-mass pre-gas- 
expulsion clusters will lose a significant fraction of their 
mass within the first s e veral M yr and rapidly expand 
(IBaumgardt fc Kroupal 120071 : iParmentier &; Gilmord 
12007ft . Since our FP model does not consider the 
effect of gas expulsion, our initial GC models are to 
be regarded as models at several Myr after cluster 
formation. 



4. SIZE EVOLUTION OF INDIVIDUAL GLOBULAR 
CLUSTERS 

The three main drivers of GC size evolution are the 
two-body relaxation, the mass loss by stellar evolution, 
and the galactic tides. In this section, we discuss the 
size evolution of individual GCs with a subset of our 
FP calculations. Figure [3] shows the ratios between Th 
values at the present time (13 Gyr) and at the beginning 
from our FP calculations as a function of r^o and Mq for 
two different Rg,o values (subscripts denote the initial 
value, hereafter). 

Two-body relaxation causes GC core to collapse 
and the subsequent formation of dynamical bina- 
ries in the core makes the whole cluster expand. 
For GCs that have undergone core collapse in the 
early phase of evolution, the size of the post-core- 
collapsc expansion follows a scaling relation oc 

M~ 1 / 3 t 2 / 3 (iGoodmanl [l98l IKim. Lee, fc Goodman! 
119981 : IBaumgardt. Hut, fc Heggid I2002I L and thus for a 
given initial mass and epoch, r/j/r^o is simply propor- 
tional to r^ . Figure [3] indeed shows that the size of 
the GCs with the same Mq tend to converge to a sin- 
gle value (rh/rh,o oc r^ ), if the GCs have small t r h,o 
(logUhfi/yr < 9). 

Mass loss by stellar evolution causes GCs to adiabati- 
cally expand to maintain virialization, and the GC sizes 
evolve following r^/r^o oc Mq/M when the stellar evo- 
lution is the main driver of the GC size evolution (jHillsl 
1980). The combination of Krou pa IMF and the stel - 
lar evolution recipe described bv iSchaller et al.l (1992) 
yields a mass loss of ~ 40% within 13 Gyr. Thus, GCs 
would expand by a factor of ~ 1.67 as a result of the stel- 
lar evolution, if two-body relaxation or the galactic tides 
are relatively less important in driving the size evolution. 
Indeed, clusters with logi^o/yr > 9 and 5ft < 0.05 have 
i"h,i3/i"h,o values between 1 and 2. 

While stellar evolution and two-body relaxation cause 
clusters to expand, galactic tides make clusters shrink in 
general. A cluster extending farther than r j (overfilling; 
r t > Tj) loses stars outside rj within a few dynamical 
timescales, and this naturally causes the mean size of the 
cluster to decrease. Since rj oc Rg(M /Mg) 1 ^ 3 where 
Mg is an enclosed mass of the Milky Way in a given 
Rg, the size decrease caused by the galactic tides takes 
place mostly while the cluster approaches R p . The clus- 
ter re-expands so mewhat by two-body relaxat ion while 
approaching R a (jBaumgardt fc Makinol I2003T ). but its 
size gradually decreases while repeating orbital motions. 
We find that clusters with 0.4 < r t /rj < 1 can also 
shrink moderately as a result of the galactic tides even 
if it underfills, and clusters initially with r t /rj < 0.4 (or 
5ft < 0.05; i.e., "isolated" GCs) can gradually move into 
the " tidal" regime as they lose mass or expand by stellar 
evolution or two-body relaxation. Figure [3] shows that 
GCs with larger 5fto are smaller at 13 Gyr for a given Mq 
and i?G,0i a s expected. 

Among various initial GC parameters, r^.o is the most 
important parameter in the size evolution caused by two- 
body relaxation (r/j/r^o oc M 1 ^ 3 r~ h 1 ) and that result- 
ing from galactic tides (5fto oc M^ 1 ^ 3 R^^rhfl for a flat 
rotation curve). For this reason, initially small GCs gen- 
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log M /M s log t rhi0 /yr 
4.0 □ 6.0 
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Fig. 3. — Ratios of values at 13 Gyr and at the beginning from some of our 1,152 Fokker— Planck calculations as a function of 
at the beginning for two different Rq,o values (4.6 kpc for the left panel and 46 kpc for the right panel) and four different initial mass 
(logMo/M© = 4, 5, 6, and 7). The approximate initial half-mass relaxation times and the initial tidal filling ratios are marked with 
different colors and symbol sizes, respectively. The blue line indicates the location of r^/r^o = 1-67, which is the expected expansion ratio 
mainly by the stellar evolution, and the red lines represents the relation r^/r^o oc F^" Q' which is the expected result when the evolution is 
dominated by the two-body relaxation. 
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Fig. 4. — Comparison of mass functions (a), radial distributions (b), and size distributions (c) at 13 Gyr (solid lines) and at the beginning 
(dashed lines) from the best-fit initial parameter set for each SD model. 



erally expand (by two-body relaxation), while initially 
large GCs generally shrink (by the galactic tides) as they 
evolve. The size evolution of intermediate GCs is deter- 
mined by more than one dynamical effect, and some GCs 
can even maintain their initial size over their whole life- 
time. 



5. SYNTHESIS OF FOKKER-PLANCK CALCULATIONS 

As discussed in Section 3, we performed a total of 1152 
FP calculations with different initial cluster conditions in 
four-dimensional parameter space, M, r^, R a , and e. The 
goal of the present study is to find the initial distribution 
of these variables that best describe the observed GCs. 



For the initial MF model, we adopt a Schechter function, 

dN(M) oc M- a exp(-M/M s )dM, (1) 

and for the initial RD model, we use a softened power-law 
function, 



dN(R G ) oc 4irR 2 G dR G /[l + {R G /R S ) 



01 



(2) 



We assume that the initial MF is independent of initial 
R G . For the sake of simplicity, we do not parameter- 
ize the distribution for e, and adopt the fixed isotropic 
distributions, i.e., dN(e) oc e de. Unlike M, the R G of 
each FP model evolves by oscillating between R p and R a . 
and thus the model RD at 13 Gyr constructed from our 
population synthesis may suffer from significant random 
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Fig. 5. — Comparison of size distributions at 13 Gyr (solid lines) 
and at the beginning (dashed lines) from SD Model 1. The upper 
panels are for initially small GCs (t r h,o < 0.5 Gyr or SRo < 0.05), 
and the lower panels are for initially large GCs (t r fe Q > 0.5 Gyr or 
SRq > 0.05). 



noise. To reduce this noise, we build a model RD by 
summing the probability distributions between R p and 
R a that are given by the orbital information at 13 Gyr, 
and we call this a phase-mixed RD. Hereafter, RDs in 
this paper refer to the phase-mixed RD. 

For initial SDs, we use six distribution models (see Ta- 
ble [T}. Models 1, 2, and 3 represent a Gaussian distri- 
bution of r/j, ph (mean density within r^), and 5R, re- 
spectively, implying that the initial GCs have the pre- 
ferred initial rh, Ph, and 5ft, with dispersions. The ini- 
tial Th of Model 1 does not correlate with the initial M 
or Rg, while Models 2 and 3 have initial correlations 

of r h (x M 1 / 3 p~ 1/3 and r h cx M^ 3 R G 2/3 ^. In Mod- 
els 4, 5, and 6, the initial is determined by pow- 
ers of initial M and/or Rg- Note that the power of 
Model 6 (r/j oc M 615 ) corresponds to that of the mass- 
size relation derived from the Faber- Jackson relation for 
early-type galaxie s (jFaber et al.lll989t iHasegan et.ll2005t 
IGieles et al.ll27)10[ ). 

Once the calculations of the 1152 FP models are done, 
the aforementioned sets of initial MF, RD, and SD mod- 
els are used to search for the best-fit parameters in five 
to seven dimensional space, depending on the SD mod- 
els (Models 1-6). For this, we synthesize our 1152 FP 
calculations with appropriate weights to produce a given 
initial MF, RD, and SD, and find a set of parameters 
that best fit the present-day MF, RD, and SD for each 
of the six SD models. When finding the best set of pa- 
rameters for each SD model, we minimize the sum of 
X 2 values from all of the L, Rg, and Rh histograms, 
which are constructed by using eight bins between 10 
and 10 5 8 Lq for L, nine bins between 10° and 10 16 kpc 
for Rg, and nine bins between 10 -0 6 and Id 1-2 pc for Rh, 
all equally spaced on a logarithmic scale. Recall that we 
use dynamical (theoretical) properties M and for set- 



ting the initial distributions, while observable quantities 
such as L and Rh are used for comparing the models and 
observations. 

6. BEST-FIT INITIAL DISTRIBUTION OF THE GALACTIC 
GLOBULAR CLUSTER SYSTEM 

The best-fit parameter sets that minimize the \ 2 values 
between observations and our calculations are presented 
in Tabled] for the six SD models. We examine the char- 
acteristics of our best initial MFs, RDs and SDs in turn. 

6.1. Initial Mass Function 

The best-fit a values for all six SD models are quite 
low, ranging between 0.01 and 0.07. The best-fit 
logM s /M Q values for all six SD models are similar to 
each other, having values between 5.8 and 5.9. Note that 
Schechter functions with such small a values are similar 
to log-normal functions, while those of a > 2 are closer 
to power-law functions. Thus, our small a values sug- 
gest that log-normal functions better describe the initial 
MF of the Galactic GC system than power-law functions 
(see Figure |4jz) , and this result is consistent with the 
result of Paper I. One way to explain the log-normal- 
like initial MF is expulsion of the remnant gas due to 
star formation in the pre-gas-expulsion cluster, which can 
quickly alter a power-law MF into a log-normal-like MF 
(jParmentier fc Gilmord 120071 ) . Another possible mecha- 
nism resulting in a rapid change in the initial MF is the 
collisions of clusters with dense clou ds or other cluste rs 
during the early phase of the galaxy (lElmegre cn 2010). 

6.2. Initial Radial Distribution 

Initial RDs from the best-fit parameter sets for all 
six SD models have similar (3 values (4.0-4.5) but a 
rather wide range of R s values (0.3-3.6 kpc), and this 
is consistent with the result of Paper I ((3 — 4.2 and 
R s = 2.9 kpc). 

Figure \4$> shows that most of the GCs that disrupt 
before 13 Gyr are located in the bulge regime (Rg,o < 
3 kpc), and most of the GCs formed in the bulge do 
not survive until now. We find that only 0.1-8.4 % of 
the total GC mass initially inside 3 kpc remains in GCs 
at 13 Gyr, and the total stellar mass that escaped from 
the GCs inside 3 kpc during the last 13 Gyr amounts to 
5 x 10 7 -3 x 10 8 M , depending on the SD model. 

6.3. Initial Size Distribution 

The initial SDs from our best-fit parameter sets are 
of larger dispersion than the present-day SDs for all six 
SD models (see Figure 3b). The initial SDs evolve into 
the narrower present-day SDs by two main effects: (1) 
expansion of GCs with small rh$, which normally have 
small t r hfl and/or small Jio, due to two-body relaxation, 
and (2) shrinkage (evaporation) of large r/j i0 GCs, which 
normally have large t r h o and/or large 5fto, due to the 
Galactic tides. Figure [Sj shows that the SDs of initially 
small GCs (upper panels) indeed shift to the larger rh 
region and those of initially large GCs (lower panels) 
shift to the smaller rh region after 13 Gyr. 

Three p- values (significance levels) for \ 2 tests of LFs, 
RDs, and SDs are acceptably high, except for Model 3, 
which has relatively small x 2 p-values for RDs and SDs 
(see Table [2]). However, the high p- values from the x 2 
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Fig. 6. — L (top left), R (top middle), and R^ (top right) histograms at 13 Gyr (thick solid lines) for SD Model 1 with the best- 
fit parameter set. Also shown together in the upper panels are the corresponding initial distributions (dashed lines) and the observed 
distributions (thin solid lines). The lower panels show the correlations between Rh and Rq (bottom left), Rh and L (bottom middle), 
and L and Rq (bottom right) relationships for the corresponding best-fit models in the upper panels (asterisks) and from the observations 
(open circles). 
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Fig. 7.— Same as Figure 5, but for SD Model 2. 
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Fig. 8. — Initial SDs of the GCs that survive until 13 Gyr in 
our best-fit SD models, Models 1 (thick solid line) and 2 (dashed 
line). Also shown is the currently observed SD (thin solid line with 
a shaded area). The overall size of the GCs were larger at birth 
than now even when only the surviving GCs are considered. 

tests do not necessarily guarantee that the models with 
the best-fit parameters restore the observed correlation 
between L, Rq, and Rh as well. Thus, we implement 
Student's f-tests to see if our models with the best-fit 
parameters agree with the observed Rq dependence of 
SDs (the Rh-Rg correlation), the L dependence of SDs 
(the Rh~L correlation) , and Rq dependence of LFs (the 
L-Rg correlation). For the Rh~RG correlation, we calcu- 
late x 2 for the difference of {log Rh) and (J\ogR h between 
the model and the observation as follows: 

X 2 «log R h ))=J2 ^ Rh '°' j) " (1 ° g Rh ' m ' j))2 



(CTl ogJ 



CTlog R h , m ,jY 



'h> S R h , mJ 



/2N , 



(3) 



where subscripts o and m stand for the observation and 
the model, respectively, subscript j represents the equal 
number Rg bins, and (. . .) denotes the averaged values. 
The same calculation is applied to Rh~L and L-Rq cor- 
relation as well. We find that Models 3-6 have f-test 
p- values that are too small (< 1%) for at least one of the 
Rfi-Rc, Rh~L, and L-Rq correlation. For this reason, 
we reject Models 3-6 as being a plausible initial SD can- 
didate. Hereafter, we call SD models 1 and 2 "the final 
best-fit SD models" . 

Figures [6] and [7] show our two remaining best-fit SD 
models, a Gaussian distribution of r/j (model 1; r/i lC = 
6.4 pc, ay h = 2.7 pc) and a Gaussian distribution of 



ph (model 2; p h . 



690 M Q pc- 



4.6 Mqpc-^). 



Note that r/^o values are not correlated with the Rg,o m 
either model. This implies that the Th,o of GCs probably 
does not depend on the strength of the galactic tides. 
Therefore, we interpret the observed, present-day Rh~Rc 
correlation (see Figure[TJ>) as an outcome of a preferential 



disruption of the larger GCs at smaller Rq due to the 
Galactic tides. 

7. DISCUSSION 

The typical Rh,o value from our final best-fit SD mod- 
els (Models 1 and 2) is ~ 4.6 pc (r^o ~ 7 pc), and 
this is 1.8 times larger than that of the present-day GCs 
(~ 2.5 pc). T his result is rathe r differ ent from a recent 
argument by iBaumgardt et al.l (|2010H that most GCs 
were born compact with r^o < 1 pc. Our result implies 
that GCs initially have a rather wide SD, the typical 
value of which is similar to that of YMCs in parsec scale, 
and have evolved to have a narrower SD with a smaller 
mean value. 

We also find that GCs formation favors a "tidal" en- 
vironment over an "isolated" environment. The number 
of tidal GCs (5R > 0.05) at Gyr from our final best-fit 
SD models is approximately five times larger than that 
of isolated GCs (9?o < 0.05). The ratio of tidal to iso- 
lated GCs, however, drastically decreases as GCs evolve 
because tidal GCs are more easily disrupted, and this 
ratio becomes ~ 0.2 at 13 Gyr. 

Figure [8] shows the initial SDs of the GCs that survive 
until 13 Gyr in Models 1 and 2. We find that these initial 
SDs are broader (a(Rh) = 2.1 and 2.5 pc, respectively) 
and centered at higher values (Rh = 4.1 and 4.0 pc) 
than the currently observed SD (a(R h ) — 1.2 pc, Rh = 
2.5 pc). Thus, the overall size of the GCs were larger at 
birth than now by a factor of ~ 2 even when only the 
surviving GCs are considered. 

The initial total masses in GCs (M t .o) of the final 
best-fit SD models are 2.8xl0 8 M (Model 1) and 
5.3xl0 8 Mq (Model 2), and the masses that have left 
the GCs during the lifetime of the Galaxy (AMt) are 
2.5xl0 8 M e (Model 1) and 5.0x10 s M (Model 2). 
These give AM t /M t>0 values of 0.89 and 0.94 for Mod- 
els 1 and 2, respectively. Our AAl t values are several 
times larger than previou s estimates made by[Baumg;ardt 
(fl998l 4.0-9.5 xlO 7 MrA IVesperinl (fl998l . 5.5 xlO 7 M ), 
and Paper I (1.5-1. 8x10 s M ). Our larger AM t values 
are due to the facts that (1) we consider virtually all dis- 
ruption mechanisms in the calculations for the dynamical 
evolution of individual GCs, and (2) we use more a flex- 
ible initial rh distribution, which can have a relatively 
larger fraction of GCs with a large 77, (larger GCs are 
more vulnerable to the galactic tide). Note that AM t 
will be larger if one considers the clusters that have been 
disrupted in the process of remnant gas expulsion. 

We note that contrary to the finding in the present pa- 
per, detailed dynamical modeling of individual clusters 
shows that at least some of the clusters must have started 
with a ver y small size. For e xamp le, M onte Carlo calcu- 
lation s by Heggie fc Gierszl (|2008h and iGiersz fc Heggid 
(|2009l 2011) find 0.58 pc, 0.40 pc, and 1.9 pc as best-fit 
initial rh values for the observed current states of M4, 
NGC 6397, and 47 Tuc, respectively. These values are 
several times smaller than the typical initial rh found for 
the Galactic GC system from our calculations, ~ 7 pc. 
However, we also note that the Monte Carlo models used 
for these three clusters all assume circular cluster orbits 
while M4 and NGC 6397 have moderate to high orbit 
eccentricities (0.82 and 0.34, respectively). We have per- 
formed several FP calculations for these two clusters and 
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find that consideration of appropriate eccentric orbits can 
increase the best-fit initial Th by a factor of 3-5. 

8. SUMMARY 

We have calculated the dynamical evolution of Galac- 
tic GCs using the most advanced and realistic FP model, 
and searched a wide parameter space for the best-htting 
initial SD, MF, and RD models that evolve into the 
present-day distribution. We found the initial MF of the 
Galactic GC system is similar to the log-normal func- 
tion rather than the power-law function, and the RD 
of the GC system undergoes significant evolution inside 
Rg = 3 kpc through the strong Galactic tides. We also 
found that the initial SD of the GC system evolves to 
narrower present-day SDs through two effects: shrink- 
age of large GCs by the galactic tides and expansion of 
small GCs by two-body relaxation. The typical initial 
projected half-mass radius from the final best-fit model, 
~ 4.6 pc, is 1.8 times larger than that of the present-day 
value, - 2.5 pc. The ratio of "tidal" GCs to "isolated" 
GCs is ~ 5 at Gyr and decreases down to ~ 0.2 at 13 
Gyr. 

Since tidal GCs are found to be dominant in the begin- 
ning, one might expect the initial size of the GCs to be 
correlated with the Jacobi radius, i.e., to be a function of 
the galactocentric radius. However, our final best-fit SD 



models (Models 1 and 2) do not seem connected to the 
galactocentric radius. This implies that the GC forma- 
tion process favors a certain size and density, regardless 
of the tidal environment. Such a i?c-independent initial 
SD evolves into a present-day SD, which shows a tight 
rh-R-G correlation through evaporation and two-body re- 
laxation. 
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SD model Functional form Parameters 

1 AT(r h>c ,<T^) a a,M a ,P,R a ,r hiC ,(TT h 

2 N{Ph,c,v 2 Ph Y a,M 3 ,l3,R 3 ,p htC ,a Ph 

3 A r (J? c ,o-|) a a,M s ,/3,i? s ,sR c ,<T S 

4 r h = KM x R v a,Ms,fS,R s ,K,X,u 

5 r h = kM x a,M a ,P,R a ,K,\ 

6 r t = nM M15 g,M s ,l3,R s ,K 



TABLE 1 
Initial SD models 



TABLE 2 

Best-fit parameters for initial GC distributions 







MF 




RD 


SD 


p-values (x 2 


test) 


p-valucs (t-test) 




Model 


a 


logM s 


P 


Rs 




LF 


RD 


SD 


Rh-Ra 


Rh-L 


L-R G 


Mt,o 






(M ) 




(kpc) 




(%) 


(%) 


(%) 


(%) 


(%) 


(%) 


(10 8 M ) 


1 


0.06 


5.8 


4.4 


1.6 


r h _ c =6.4, <r rh =2.7 


72 


37 


87 


31 


92 


48 


2.8 


2 


0.01 


5.9 


4.2 


0.3 


p hiC =690, o- ph =4.6 


66 


43 


63 


28 


42 


22 


5.3 


3 


0.01 


5.8 


4.5 


3.6 


5R C =0.04, (7^=0.02 


94 


8 


11 


1 


79 


38 


1.1 


4 


0.01 


5.8 


4.0 


1.9 


K =10- 31 , A=0.45, 1^=0.3 


86 


27 


43 


3 


1 


56 


1.3 


5 


0.07 


5.9 


4.0 


1.7 


K =10- 2 - , A=0.44 


24 


28 


25 


1 


1 


10 


1.9 


6 


0.05 


5.9 


4.4 


1.9 


K =10- 2 - 9 


74 


21 


68 


15 








5.4 



Note. — The p value for the x 2 test (t-test) is the probability of having a x 2 (t) value that is larger than the value obtained from our x 2 
(t) test between the model and the observation, whose degree of freedom is 8 or 9 (4). rh and rh, c are in units of pc, and ph and p^.c are in 
units of M /pc 3 . 



